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ABSTRACT 

Complete  solutions  to  the  equations  of  motion  are  possible  using 
numerical  schemes.  The  model  developed  by  Gait  (1972)  was  investigated 
for  stability  and  convergence.  A  series  of  runs  indicated  that  the 
scheme,  the  three  level  Adams-Bashforth  method  for  the  new  value  of 
vorticity  and  then  the  successive  over-relaxation  method  for  the  stream 
function,  was  stable  for  transport  not  exceeding  100  Sverdrups  and  a 
time  step  of  one-half  pendulum  day.  Convergence  was  determined  by 
comparisons,  with  known  analytic  solutions,  of  phase  velocity  of  Rossby 
waves,  and  for  the  steady  state  conditions,  patterns  of  the  stream 
function  and  the  vorticity.  Improved  accuracy  of  the  model  was  attained 
by  using  a  finer  grid. 
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I.   INTRODUCTION 

The  study  of  the  dynamics  of  oceanic  flow  has  been  attempted  by 
using  the  equations  of  motion.  Since  solutions  to  those  complex  equa- 
tions in  their  complete  form  was  not  probable  it  became  necessary  to 
simplify  them  by  making  highly  restrictive  assumptions  that  allowed 
for  some  of  the  functional  relationships  to  be  defined  as  constants  and 
the  non-linear  terms  to  be  disregarded.  Ekman  (1905)  developed  the  now 
classic  boundary  layer  solution,  known  as  the  Ekman  spiral,  which 
showed  the  effects  of  wind  stress  acting  on  a  free  surface.  Stommel 

(1948)  was  able  to  reduce  the  complex  equations  to  a  balance  between 
wind  stress,  a  parameterization  of  friction,  and  Coriolis  force.  His 
solution  dramatized  the  phenomenon  of  western  intensification. 

The  derivation  of  the  vorticity  equation  from  the  vertically  inte- 
grated equations  of  motion  removed  the  necessity  of  defining  the  press- 
ure field  at  the  surface  and  simplified  the  boundary  conditions.  Munk 

(1949)  developed  a  solution  to  a  simplified  form  of  the  vorticity  equa- 
tion which  contained  no  non-linear  terms,  was  invariant  with  time, 
assumed  a  zonal  wind  stress,  and  assumed  a  constant  eddy  coefficient. 
In  both  solutions  the  features  Stommel  observed  in  his  solution  were 
evidenced  and,  additionally,  Munk  (1949)  observed  the  counter  current 
region  as  well . 

Development  of  analytic  solutions  to  the  equations  of  motion 
necessarily  demanded  regular,  closed  boundaries,  and  for  some  forms,  a 
constant  depth.  In  each  solution  mentioned,  an  ocean  of  constant  depth 
was  assumed  or,  as  in  Ekman's  solution,  the  ocean  was  considered  as 


unbounded  in  depth.  To  study,  analytically,  ocean  basins  of  irregular 
boundaries,  with  inflow  and  outflow,  and  variable  depth  lead  to  an 
impossible  task. 

Early  work  in  developing  these  solutions  was  usually  attempted  for 
ocean  basins  that  were  rectangular  in  shape.  A  need  to  investigate  the 
effects  of  a  laterally  sloping  boundary  was  recognized,  as  evidenced 
by  the  more  recent  work,  by  Munk  and  Carrier  (1950)  which  was  done 

on  a  triangular  shaped  basin.  Studies  of  the  effects  of  bottom  topo- 
graphy on  ocean  circulation  have  been  presented  by  Warren  (1963)  and, 
more  recently,  by  Veronis  (1966).  The  solutions  Veronis  developed 
were  to  the  time  dependent,  non-advective  equations  with  no  friction  or 
wind  stress  and  fora  rectangular  basin. 

It  was  clearly  indicated  that  with  each  solution  a  different  pheno- 
menon of  oceanic  flow  could  be  described  using  the  vorticity  equation 
or  the  equations  of  motion.  The  optimum  stiuation,  then,  would  be  to 
obtain  a  complete  solution  to  these  equations  so  that  instead  of  being 
restricted  to  solutions  which  describe  only  one  of  the  phenomenon  a 
complete  interaction  of  advection,  Coriolis  force,  friction,  and  wind 
stress  could  be  studied  as  a  function  of  time. 

Since  a  numerical  solution  to  the  equations  of  motion  removed  the 
necessity  of  simplification  and  allowed  for  a  complete  solution,  various 
models  have  been  developed.  Usually  the  only  simplifications  involved 
were  to  specify  the  density  field  and  friction  as  constant.  More 
recently  Gait  (1972)  described  development  of  a  numerical  model  of  a 
constant  density  form  of  the  vorticity  equation. 


The  non-dimensional   equations  that  were  used  are  represented  by 
equations   (1)  and  (2). 
Vorticity: 

M+h7.,H(s*i).«,get  vff) 


(i) 


Stream  function: 

1 
h 


V.(rVM    =   5  (2) 


Where: 

5  =  vorticity 

h  =  depth 

f  =  Coriolis  parameter 

K  =  lateral  kinematic  eddy  coefficient 

~  =  wind  stress 

vH  =  horizontal   gradient  operator  stream  function 

In  this  form  the  non-dimensional ized  constants  a   and  6  were  defined: 

a  =     ^Q/FL2D 
6   =     K/FL2 

Where: 

V  =  a  representative  value  of  the  stream  function 

F  =  a  representative  value  of  the  Coriolis  parameter 

D  =  mean  depth 

L  =  the  grid  length 

The  non-dimensional  wind  stress  t  was  related  to  x',  the  dimensional 

form,  by  equation  (3). 

/  p  F  ^  \ 
*  =  T'(— TT^j  (3) 


Where: 

p    =      density 
The  non-dimensional   stream  function  was  then  defined: 

*     -     *'    (FL2D)     =     *'    ♦ 

o 

Using  the  non-dimensional   forms  of  the  equations  the  model  was 
constructed.     The  new  value  of  vorticity  was  obtained  by  the  three  level 
Adams-Bashforth  differencing  method. 

The  stream  function  was  then  relaxed  to  within  acceptable  limits  by 
the  method  of  successive  over-relation.     This  scheme  was  applied  to  a 
polygonal   computational  molecule. 

This  computational   grid  allowed  for  the  specification  of  irregular 

boundaries  and  for  inflow  and  outflow.     Newtonian,  lateral   friction  was 

2 

parameterized  in  terms  of  the  horizontal  Ekman  number,  K/FL  . 

For  models  which  have  been  developed  it  has  been  necessary  to,  first, 
determine  their  stability  and  second,  to  determine  their  convergence  to 
the  true  solution.  Since  complete,  analytic  solutions  were  not  available, 
these  models  were  run  for  conditions  which  were  representative  of 
simplified  equations  to  which  solutions  existed.  Comparisons  were  then 
made  to  assess  the  validity  of  the  models. 

For  some  finite  difference  schemes  it  is  possible  to  mathematically 
solve  for  errors  generated  and  determine  rates  of  convergence  as  well  as 
showing  what  bounds  there  must  be  on  the  spatial  and  temporal  steps  to 
ensure  stability.  In  referring  to  the  model  Gait  developed  it  was  not 
possible  to  make  general  calculations  for  convergence  and  stability. 
As  indicated  by  Baer  and  Simons  (1970)  a  method  to  investigate  numerical 
stability  in  a  model  can  be  effected  by  monitoring  a  conservative 
property  such  as  energy  and  total  vorticity.  The  method  used  to  determine 

10 


convergence  was  by  a  comparison  of  the  patterns  of  vorticity  and  the 
stream  function  with  known  analytical  solutions.  This  method,  for 
example,  was  used  by  Blandford  (1970). 

Establishment  of  the  characteristics  of  Gait's  model  in  terms  of 
stability  and  convergence  was  done  in  two  parts;  1):  For  stability,  a 
form  of  kinetic  energy,  total  vorticity,  and  the  maximum  value  of  the 
stream  function  were  monitored;  2):  To  check  convergence  a  comparison 
was  made  between  known  analytic  solutions  and  the  solutions  predicted 
by  the  model.  These  comparisons  used  patterns  of  the  vorticity  and  the 
stream  function,  and  of  the  phase  velocity  of  Rossby  waves.  These 
comparisons  were  made  with  two  different  model  configurations. 

Once  the  basic  shape  and  boundary  conditions  were  determined  the 
initial  input  was  either  constant  or  variable  depth,  distribution  of 
the  Coriolis  parameter,  magnitude  and  direction  of  the  wind  stress, 
stream  function  and  vorticity  patterns,  and  the  local  rate  of  change  and 
vorticity  one  time  step  back.  The  parametric  specification  of  density 
and  frictional  coefficients,  which  were  assumed  constant,  and  the 

physical  dimensions  of  a  particular  configuration  were  determined  by 

?         2 

the  magnitudes  of  the  non-dimensional izing  constants,  i>  /FL  D  and  K/FL  . 

Additionally,  for  those  runs  when  these  constants  were  zero  the 
dimensions  of  the  model  were  determined  by  the  distribution  of  the 
Coliolis  parameter,  magnitude  of  the  wind  stress,  and  bottom  topography. 

The  first  configuration  was  a  rectangle  shape  that  could  only  be 
approximated  by  the  model  due  to  the  use  of  the  polygonal  molecules 
for  the  computational  grid.  The  boundary  conditions  were  set  up  so  as 
to  represent  an  enclosed  basin  with  free  slip  boundaries  with  no 
friction  in  the  model.  The  non-linear  terms  were  not  used  and  the  wind 
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stress  set  to  zero  so  that  this  solution  could  be  compared  to  the  free 
wave  solutions  of  Veronis  (1966).  One  run  was  with  constant  depth  and 
a  variable  Coriolis  parameter  and  the  other  run  was  with  variable 
depth  and  a  constant  Coriolis  parameter. 

The  second  series  of  runs  were  done  using  a  parallelogram  shape 
which  had  closed  regular  boundaries.  The  east-west  boundaries  were 
inclined  sixty  degrees  from  the  horizontal,  measured  in  a  positive 
sense.  This  was  run  with  constant  depth  and  a  variable  Coriolis 
parameter.  The  north-south  boundaries  were  set  to  free  slip  conditions 
and  the  east-west  boundaries  were  set  to  no  slip  conditions  with 
friction  in  the  interior.  A  zonal  wind  stress  set  in  as  a  consine  func- 
tion with  maximum  values  at  the  north-south  boundaries.  The  model  did 
not  accurately  convert  this  to  the  curl  of  wind  stress  and,  therefore, 
all  successive  runs  were  made  with  a  direct  input  of  the  curl  of  the 
wind  stress.  These  particular  requirements  were  so  that  a  comparison 
with  a  modified  solution  of  Munk  and  Carrier's  (1950)  could  be  made. 

The  investigation  of  the  effect  of  reducing  the  spatial  step  size 
was  done  for  both  configurations  by  increasing  the  number  of  points 
and  appropriately  scaling  the  input  parameters.  Reduction  of  the  grid 
length  to  one-half  was  accomplished  so  as  to  better  define  the 
velocities  and  western  boundary  layer. 
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II.  STABILITY  ANALYSIS 

Baer  and  Simons  (1970)  indicated  that  for  the  Adams-Bashforth 
iteration  scheme,  when  used  on  a  non-linear  equation,  it  was  sufficient 
to  determine  stability  for  a  linear  form  of  that  equation  since  the  on- 
set of  instability  appears  to  be  the  breakdown  of  the  linear  solutions. 
Since  determination  of  computational  stability  was  desired  it  was  suffi- 
cient, (Baer  and  Simons  [1970]),  to  monitor  the  conservative  properties 
which  were,  kinetic  energy,  total  vorticity,  and  the  maximum  value  of 
the  stream  function.  For  the  stable  case,  Figure  1,  these  properties 
were  bounded  over  the  period  of  integration.  Figure  2,  the  unstable 
case,  indicates  some  of  these  values  as  growing  unbounded  after  a  certain 
integration  time.  The  model's  response  for  the  stable  situation  was  a 
comparatively  short  spin-up  time,  as  indicated  by  the  values  of  the 
properties  monitored.  The  kinetic  energy  and  maximum  value  of  the 
stream  function  reached  a  maximum  value  at  this  time  and  were  then 
observed  to  oscillate  throughout  the  remainder  of  the  integration  time. 
This  response  was  similar  to  the  results  reported  by  Blandford  (1970) 
for  a  model  in  which  was  specified  a  lateral  viscosity  and  no  slip 
boundary  conditions.  Total  vorticity  reached  as  maximum  value  later  in 
the  integration  then  remained  constant  throughout  the  remainder  of  the 
integration  time.  For  the  unstable  situation  the  kinetic  energy  after 
a  period  of  time  began  to  grow  in  an  exponential  manner  with  oscilla- 
tions of  the  stream  function  increasing.  The  vorticity,  for  this 
situation,  grew  steadily  throughout  the  integration  time. 
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To  more  accurately  determine  the  limits  of  stability  various 
conditions  were  imposed  on  a  constant  depth  model  basin  with  basic 
configuration  as  shown  in  Figure  3. 


Wind  Stress 


N 


->x 


Figure  3.  Basic  Model  Configuration  for  Stability  Analysis 

For  all  but  two  runs  mixed  boundary  conditions  were  used  with  the  north 
and  south  boundaries  being  free  slip.  Of  the  other  two  runs,  one  was 
run  with  free  slip  boundary  conditions,  the  other  with  no  slip  boundary 
conditions. 

The  distribution  of  the  wind  stress  over  the  model  basin  was  zonal 
in  nature.  The  direction  of  the  wind  stress  was  the  cosine  function  in 
the  north-south  direction  with  the  maximum  stress  in  the  westward 
direction  at  the  south  boundary  and  the  maximum  eastward  stress  at  the 
north  boundary. 

The  first  series  of  runs  were  made  with  the  constant  (^0/FL  D)  set 
to  zero  which  resulted  in  the  linear  form  of  the  equations.  The  wind 
stress  was  set  at  a  representative  value  and  the  Coriolis  parameter  was 
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allowed  to  vary  over  the  constant  depth  model  basin.  The  physical  dim- 
ensions of  the  basin  were  determined  from  this  distribution.  Table  I 
shows  the  stability  for  a  range  of  friction,  K.  For  the  grid  length 

L  =  4.6  x  10  m,  F  =  10"  sec"  ;  K  ranges  from  approximately  7.5  x  10  to 

7  2-1 

7.5  x  10  cm  sec  .  The  effect  of  varying  the  magnitude  of  the  wind  stress 

for  this  test  had  the  effect  of  varying  the  transport,  and  by  using 
large  values  for  the  wind  stress  the  model  was  unstable. 

Table  I.  Stability  Dependence  on  Friction  and  Time  Step. 


K/FL2 

DT 

KDT/FL2 

STABLE 

2.86  x 

io-3 

0.5 

1.43  x  10"3 

Yes 

2.86  x 

10"4 

0.5 

1.43  x  10"5 

Yes 

2.86  x 

10"4 

1.0 

2.86  x  10"4 

Yes 

2.86  x 

10"4 

2.0 

5.72  x  10"4 

Yes 

2.86  x 

10"4 

3.0 

8.58  x  10"4 

? 

2.86  x 

10"4 

5.0 

1.43  x  10"3 

? 

2.86  x 

10"4 

7.5 

1.71  x  10'3 

No 

2.86  x 

10"4 

10.0 

2.86  x  10"3 

No 

8  ?    1 
Allowing  K  to  remain  constant  at  7.5  x  10  cm  sec"  the  wind  was 

varied  as  shown  in  Table  II. 

Table  II.  Stability  Dependence  on  Transport 


curl  t  DT  STABLE 

-sin(y)  0.5  No 

0.0  0.5  Yes 

-10"4sin(y)  0.5  Yes 

-10"2sin(y)  0.5  Yes 

-10"]sin(y)  0.5  Yes 
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By  using  the  non-dimensional izing  coefficient  for  the  wind  stress, 

IP   F  f. 


Tx  =  Tx  [— —J 


where  p  is  density,  F  a  representative  value  of  the  Coriolis  parameter, 

ft   ^ 
L  the  grid  length,  i>     in  Sverdrups  (10  m  /sec)  was  determined.  The 

maximum  absolute  value  of  the  non-dimensional  wind  stress,  t,  was  1.0. 

with  p  =  10   grams  m"  ,  F  =  10  sec"  ,  L  =  4.6  x  10  m  the  runs  depicted 

5 
in  Table  II  represent   a  range  of  maximum  transport  from  zero  to  10 

Sverdrups. 

The  effects  of  non-linearities  on  the  stability  of  the  model  are 

7      2 

shown  in  Table  III.  Keeping  K  constant  at  7.5  x  10  ,  <P/FL  D  was  allowed 

to  vary  with  varying  wind  stress. 

Table  III.  Stability  Dependence  on  Non-linear  Terms  and 
Time  Step. 


curl  t 

^0/FL2D 
10"2 

DT 
1.0 

STABLE 

-10_1sin(y) 

Yes 

0.0 

10"2 

1.0 

Yes 

-lO^sinCy) 

10"2 

1.0 

No 

-lO^sinty) 

10"3 

1.0 

Yes 

-10sin(y) 

10"3 

5.0 

• 

-10_1sin(y) 

10"6 

1.0 

Yes 

-10"]sin(y) 

10"3 

10.0 

No 

2 
Inclusion  of  the  non-linear  terms  by  ^0/FL  D  where  D  is  the  mean  model 

-2     2 
depth,  \\>     ranged  from  10   to  10  Sverdrups.  This  also  scaled  the 

-13 

magnitude  of  the  wind  stress  and  it  ranged  from  approximately  2  x  10 

-9  2   -2 

to  2  x  10   m  sec  for  its  maximum  values.  Similar,  unstable,  response 

for  large  transport  was  obtained  by  setting  i>   /FL  D  to  10.0  which 

5 
represented  approximately  10  Sverdrups. 
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Curl 

K/FL2 

-10"4sin(y) 
-10"4sin(y) 

5  x  10"4 
5  x  TO"4 

-10"4sin(y) 

5  x  10"4 

Table  IV  indicates  the  final  three  runs  made  for  determination  of 

computational  stability.  Again,  using  the  linear  form  of  the  equations 

o  o   -1 
and  K  approximately  10  cm  sec  ,  the  maximum  transport  for  the  wind 

stress  used  was  approximately  10  Sverdrups. 

Table  IV.  Stability  Dependence  on  Time  Step. 

DT  STABLE 

1.0  Yes 

5.0  ? 

10.0  No 

An  attempt  to  compute  computational  stability  was  made  by  using; 

(c  &    +  K  -^~j)<     1.0 
\    AX        (AX)2/ 

The  term,  C  —  where  C  is  a  characteristic  velocity,  was  not  important 

AX 

in  the  series  of  linear  runs.  The  term  K  /A  \^  was  computed  as  shown 
in  Table  I  with  the  results  of  a  somewhat  larger  bounds  for  At.  For 

the  transport  range  of  10  to  100  Sverdrups  and  a  coefficient  of  friction 

8  2-1  4 

on  the  order  of  10  cm  sec   a  time  step  of  10  seconds  ensured  stability. 

This  represented  a  one-half  pendulum  day.  The  effects  of  the  different 

boundary  conditions  mentioned  previously  was  to  change  only  the  model 

spin  up  time.  This  observation  was  similar  to  that  reported  by  Blandford 

(1970). 
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III.  ROSSBY  WAVE  PROPAGATION 

Convergence  of  the  model  to  a  steady  state  analytic  solution  was 
determined  by  allowing  the  model  to  generate  time  dependent,  free  wave 
solutions  and  allowing  the  transients  to  die  out  leaving  only  the  steady 
solution.  These  Rossby  waves  can  occur  in  an  ocean  basin  and  are  of 
two  types.  The  classical  Rossby  wave  can  exist  in  a  constant  depth 
ocean  due  to  the  change  in  the  Coriolis  parameter  with  latitude.  Con- 
sidering the  situation  where  there  is  no  friction  and  no  wind  stress, 
with  the  Coriolis  parameter  increasing  from  south  to  north,  a  Rossby 
wave  will  tend  to  propagate  along  lines  of  a  constant  Coriolis  parameter. 
For  this  case,  then,  the  wave  would  propagate  in  an  east  to  west  direction 

The  second  type  of  Rossby  wave  that  can  exist  in  an  ocean  basin  is 
the  topographic  Rossby  wave.  Consider,  now,  a  constant  Coriolis  para- 
meter but  with  a  simply  varying  depth  which  increases  from  east  to 
west.   The  topographic  Rossby  wave  will  propagate  as  before  except  now 
along  lines  of  constant  depth  in  an  east-west  direction.  For  a  simply 
varying  depth,  increasing  from  east  to  west  the  topographic  Rossby 
waves  will  propagate  from  south  to  north  across  the  ocean  basin.  These 
waves  can  be  described  by  solutions  to  the  potential  vorticity  equation, 


|i  +  V-V 


(£j=;  c  =  vorticity  (4) 


Considering  that  potential  vorticity  is  conserved  following  a  parti- 
cle of  water',  equation  (4)  states  that  for  free  wave  solutions  the  second 
term  is  zero,  tiat  is,  f/h  is  constant.  Using  equation  (4),  Veronis 
(1966)  developed  a  time  dependent,  harmonic,  constant  depth  solution 
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for  the  stream  function,  y   (x,y,t). 

i>  =  \p     sin(  Ax  +  cot)  sinax  sinyy  (Veronis[1966])       (5) 

Equation  (5)  describes  sinusoidal  waves  propagating  across  a  basin  in 
the  -x  direction  enclosed  in  a  sine  envelope.  For  the  first  mode  in  the 
x  and  y  directions  the  following  constants  were  defined: 

A  =    3/2 

a  =     tt/L 

y  =     tt/M         , 


w 


3/2[a2  +  Y2]2 


gf 

3  =    constant  =  — 

L  =    basin  dimension  in  x  direction 

M  =    basin  dimension  in  y  direction 

For  use  in  initializing  the  model  non-dimensional ization  of  equation  (5) 
was  necessary.  By  non-dimensional izing  L,  M  and  x  with  a  grid  length, 
aL,  and  oo  with  F,  the  following  non-dimensional  constants  were  defined: 

X   ■     [a  l   +Y   I7 
a'   =     (tt/L)  AL 

Y'   =     (tt/0.866L)  AL 

i  i 

w   =     3AL/2X  F 

AL  =    grid  length 

L   =    basin  dimension 

F   =    a  representative  value  for  the  Coriolis  parameter 

The  phase  velocity  for  this  solution  is  co/A  which  means  that  the  non- 
dimensional  phase  velocity  is  w'A'  =  3AL/2x2  F,  or  SAL/2F[a  2  +  y  ]. 
•  A  second  solution  developed  by  Veronis  (1966)  was  for  a  linearized 
form  of  equation  (4),  but  with  a  constant  Coriolis  parameter  and  with 
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-kx 
variable  depth  that  was  assumed  to  be  h  =  he   .  This  solution, 

equation  (6),  described  sinusoidal  wave  patterns: 

$  =  $  e  '  sin(Ay  +  wt)  sinax  sinyy  (6) 

propagating  in  the  positive  y  direction  which  are  enclosed  in  a  sine 
envelope.  The  envelope  is  increasing  in  an  exponential  manner  in  the 
x  direction.  Again,  for  the  first  mode,  the  dimensional  parameters  were 
defined. 

X   =    fk/2w 

a   =     tt/L 

y   =     tt/M  1 

a)   =     fk/(4[a2  +  y2]  +  k2)2 

f   =    constant  Coriolis  parameter 

L   =    basin  dimension  in  x  direction 

M   =    basin  dimension  in  y  direction 
Using  Al_  to  non-dimensional ize  L,  M,  and  A,  and  w  with  F  the  following 
parameters  were  defined: 

k'  =     kAL  1 

i  i  o    i  o       '9  9" 

u   =      k'/(4[a  C   +Y  *]  +  k  ')' 

a'   =     (tt/L)  AL 

U/0.866L)  AL 


i 


The  non-dimensional  phase  velocity, u  /A  ,  was  obtained: 

1 

ii         i       io     19       '9  9" 

0)  /A   =  2k  /(4[a  ^  +Y  ]  +  k  d)d 

The  rectangle  shaped  basin  which  was  desired  could  not  be  exactly 
specified  as  depicted  in  Figure  4.  This  model  basin  shape  had  sawtoothed 
east-west  boundaries  due  to  the  computational  grid  being  made  up  of 
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polygonal  molecules.  Scaling  the  model  basin  to  physical  dimensions  was 
done  in  the  following  manner.  For  the  constant  depth  configuration  the 
model  was  scaled  on  the  distribution  of  the  Coriolis  parameter. 


AM 


ft  =  IE  =  AL_     n  =  9  19 

3   3y   nAM  »   n   y>  ,y     ,    , 

AF  =  0.75  x  10"Hsec"' 
a  =  60° 


using  a  representative  value  for  3  = 


1.8  x  10"11  "Isec"1, 


AM  = 


AN  = 


.75 


(9M1.8X10-11) 


=  4.63  x  10°  meters 


AM 


sin  60' 


=  5.35  x  10  meters 


The  configuration  for  which  the  Coriolis  parameter  was  constant  and 

the  depth  varied  as  in  equation  (6)  determined  the  physical  dimensions 

of  the  model  in  the  following  manner: 

h  =  h  e'kx 
o 


In 


x  = 


=  n  A  N 


n  =  9,19 

k  =  2.2  x  lO'V1 

hA  =  3.0  x  103m 

0         3 

h  =  1.0  x  10Jm 


AN  =  5.56  x  10  m 


AM  =  sin  60°(aN)  =  4.83  x  10  rn 

Since  free  wave  solutions  were  desired,  all  of  the  boundary  condi- 

2        ? 
tions  were  set  at  free  slip.  The  constants  K/FL  and  i>   /FL  D  were  set 

to  zero  so  that  the  model  had  no  friction  and  so  the  linear  equations 

could  be  solved.  Since  the  wind  stress  was  also  zero  for  these  runs 

initialization  of  the  stream  function  and  corticity  patterns  was 
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necessary.  To  accomplish  this,  the  non-dimensional i zed  solutions 
developed  earlier  from  equations  (5)  and  (6)  were  solved  at  the  non- 
dimensional  time  5  and  the  resulting  patterns  used  as  input  to  initialize 
the  model . 

In  general ,  the  model  propagated  the  two  types  of  Rossby  waves  across 
the  basin  with  the  kinetic  energy  and  total  vorticity  remaining  essen- 
tially constant.  The  fluctuations  in  kinetic  energy  that  were  observed 
had  their  maximum  value  when  the  wave,  depicted  in  Figure  4,  reached 
the  center  of  the  basin.  This  was  expected  since  the  waves  were  pro- 
pagating in  a  sine  envelope  which  reached  its  maximum  at  the  center  of 
the  basin.  The  influence  of  the  boundaries  on  the  waves  appeared  to 
cause  them  to  disperse  laterally.  This  observed  influence  quickly 
disappeared  as  the  waves  moved  to  the  interior  of  the  model 


Direction  of  Movement 


Figure  4.  Example  of  Rossby  Wave  Propagation 
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Reducing  the  grid  length  to  one-half  its  former  value  by  doubling  the 
number  of  points  on  a  side  did  not  change  these  observed  characteristics 
Subsequent  analysis  indicated  that  these  waves  were  more  accurately 
defined  with  the  decreased  grid  length. 

The  phase  velocity  for  both  series  of  runs  was  determined  by  using 
Figures  5  through  8.  The  first  comparisons  with  analytic  solutions  are 
shown  in  Table  V.  This  was  for  the  configuration  denoting  a  basin  of 
constant  depth  and  with  variable  Coriolis  parameter. 

Table  V.  Comparison  of  Rossby  Wave  Velocity. 


NUMBER  OF  GRID 
POINTS/SIDE 

GRID  LENGTHS/10 
TIME  STEPS 

NON-DIMENSIONAL  PHASE  Vel 
Model  Analytic  Solution 

10 
20 

1.2 
2.4 

.12         .169 
.24         .356 

The  percent  error  for  ten  grid  points  was  29  %   and  that  for  twenty  grid 
points  was  17  %. 

The  configuration  of  a  basin  with  an  exponentially  varying  depth 
and  with  a  constant  Coriolis  parameter  was  analyzed  in  a  similar  manner 
for  topographic  Rossby  waves.  Table  VI  indicates  the  comparison  made 
with  the  analytic  solution. 

Table  VI.  Comparison  of  Topographic  Rossby  Wave  Velocity. 

NUMBER  OF  GRID     GRID  LENGTHS/10      NON-DIMENSIONAL  PHASE  Vel 
POINTS/SIDE        TIME  STEPS        Model   Analytic  Solution 

10  1.2  .12  .212 

20  3.5  .35  .445 
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The  same  trend  was  observed  in  that  the  percent  error  for  ten  grid 
points  was  33  %  and  that  for  twenty  grid  points  was  21  %. 

This  series  of  runs  established  a  very   important  trend  for  the 
model.  Increasing  the  number  of  computational  points  increases  the 
accuracy  of  the  model  for  the  free  wave  solutions.  In  approximate 
terms,  doubling  the  number  of  computational  points  very   nearly  halves 
the  percent  error  as  compared  to  the  known  analytic  solutions.  This 
determination  was  subject  to  graphical  interpretation  which  introduced 
possible  errors.  In  an  attempt  to  clarify  the  picture,  the  figures 
were  normalized  with  respect  to  the  sinax  envelope  but  no  significant 
improvement  resulted. 
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IV.  THE  STEADY  STATE  SOLUTION 

The  assumption  that  steady  state  conditions  exist  in  an  ocean  basin 
was  made  by  Munk  and  Carrier  (1950)  when  they  derived  a  beta  plane  solu- 
tion on  a  triangular  basin  of  constant  depth.  This  solution  to  the 
vorticity  equation  balanced  frictional  effects,  Coriolis  force,  and  an 
assumed  wind  stress  distribution.  To  make  this  compatable  with  the 
configuration  of  the  numerical  model,  Figure  9,  the  dependence  of  the 
basin  width  in  the  x  direction  on  y  was  removed  so  as  to  make  it  a 
constant  width.  The  solution  was  then  corrected  for  the  boundary 
conditions  which  were  no  flow  normal  to  the  boundaries.  This  solution, 


Wind  Stress 


Figure  9.  Basic  Model  Configuration 


in  non-dimensional  form  is  shown  in  equation  (4); 

i>   =  siny[-1.158{(l-e)cos(^pa  +  .524)  -esin(^|  pa)}e("pa/2) 

+  1  -  e(pa  -  e(Pa"1/e)  +  1)]  (7) 
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with  the  following  parameters  defined 


e 

= 

y/P 

P 

'  = 

(l-*2  tan2e)'2/3 

e 

= 

60° 

$ 

= 

1-1/3(4*-  <j>0)  tan* 

yo 

= 

latitude  at  which  y  = 

0 

<J> 

= 

latitude  at  which  y  = 

IT 

a 

= 

x/y  -  (y/Y)  tane 

ya 

B 

width  of  the  basin 

The  modified  solution  is  shown  in  Figures  10  and  11.  These  patterns 
were  then  compared  to  the  steady  state  reached  by  the  model .  To 
parallel  the  model  with  this  solution  the  boundary  conditions  were 

set  at  free  slip  conditions  for  the  north  and  south  boundaries  and  no 

2 
slip  for  the  east  and  west  boundaries.  The  constant,  K/FL  ,  curl  of 

the  wind  stress,  and  the  distribution  of  the  Coriolis  parameter  were 

scaled  with  respect  to  each  other  so  as  to  represent  a  constant  depth 

3 
basin  with  approximate  physical  dimensions  of  5  x  10  km  in  the  x 

3 
direction  and  4.5  x  10  km  in  the  y  direction,  with  the  curl  of  the  wind 

-4 
stress  being  distributed  as  10  sin(y)  over  the  basin.  The  representa- 

8  2 
tive  value  for  friction,  K,  using  these  values,  was  10  cm  /sec. 

As  indicated  by  Figure  12  this  configuration  reached  steady  state 

at  about  1  x  10  seconds.  This  time  is  when  the  maximum  value  of  the 

stream  function  was  essentially  at  its  steady  state  value  with  the  vor- 

ticity  and  kinetic  energy  within  95  %   of  their  steady  state  values.  To 

determine  the  convergence  of  this  to  the  analytic  solution,  patterns  of  the 

vorticity  and  of  the  stream  function  were  drawn,  Figures  13  and  14.  The 

difference  between  these  steady  state  patterns  that  the  model  predicted 
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and  that  of  Munk  and  Carrier  could  probably  be  attributed  to  the  differ- 
ence between  the  physical  conditions  which  brought  the  two  solutions  to 
a  steady  state.  A  numerical  model  which  Gates  (1970)  used  to  solve  a 
slightly  different  form  of  the  equations  of  motion  produced  similar  steady 
state  patterns.  He  attributed  this  to  Rossby  wave  reflections  from  the 
boundaries  with  the  most  important  reflection  being  the  one  from  the 
waves  which  propagated  from  east  to  west.  These  waves  diminished  in 
magnitude  with  each  crossing  as  well  as  diminishing  in  magnitude  as  they 
propagated  through  the  model.  The  basic  differences,  not  considering 
the  finite  difference  methods,  between  Gates'  model  and  Gait's  model 
were:  all  boundaries  no  slip  as  opposed  to  mixed  free  slip/no  slip 
boundaries,  and  free  surface  as  opposed  to  a  rigid  lid.  These  differ- 
ences did  not  cause  the  two  models  to  vary  to  any  great  extent  other 
than  a  longer  spin  up  time  and  a  phase  shift  in  Gait's  model.  The  most 
important  parallel  is  that  of  the  type  of  friction  which  was  lateral 
and  did  not  include  any  bottom  friction. 

Blandford  (1970),  after  considering  different  combinations  of 
boundary  conditions,  was  able  to  show  by  using  either  free  slip  or  no 
slip  lateral  boundaries  that  bottom  friction,  expressed  as  an  exponen- 
tial boundary  layer,  dissipated  more  energy  than  did  lateral  friction. 
By  including  this  boundary  layer  in  the  solution  Blandford  noted  that 
the  Rossby  waves  which  were  generated  with  model  spin  up  tended  to  be 
damped  out  so  that  they  did  not  dominate  the  interior  flow.  The  boundary 
layer  also  caused  a  reduction  in  the  time  at  which  the  model  reached 
steady  state  conditions. 

The  response  of  this  model,  early  in  the  integration,  was  to  develop 
maximum  transport  near  the  center  of  the  basin.  This  cell  soon  propagated 
to  the  western  edge  where  it  dominated  the  flow  as  the  boundary  layer. 
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After  the  boundary  layer  was  established,  Rossby  waves  began  to  initial- 
ize in  the  southeast  corner  and  propagate  to  near  the  western  boundary. 
These  inertia!  waves  remained  near  the  southern  boundary  and  became 
established  in  the  flow  as  the  weaker  counter  current  region.  Gates 
reported  observations  on  these  waves  obtaining  the  phase  velocity  of 
5.30  meters/sec  for  their  westward  component.  These  waves,  depicted  in 
Figures  15  and  16,  had  as  their  western  component  of  the  phase  velocity 
5.34  meters/sec  for  the  larger  grid  and  4.83  meters/sec  for  the  smaller 
grid. 

To  better  define  the  phenomenon  observed  in  Gait's  model,  a  reduc- 
tion of  the  spatial  step  size  was  attempted.  Already  determined,  in 
Section  II,  was  the  improved  accuracy  of  the  solution  with  a  smaller 
grid  length.  The  response  of  the  model  for  the  stream  function  for 
different  grid  lengths,  AL,  was  observed  in  Figure  17  as  being  a  higher 
maximum  value  for  (1/2) AL.  This  difference  was  not  observed  in  the 
comparison  of  kinetic  energy,  Figure  18.  The  difference,  in  the  stream 
function,  between  1 AL  and  (1/2) AL  was  attributed  to  the  influence  of 
transients  for  the  larger  grid  length.  A  profile  of  the  stream  func- 
tion for  the  two  grid  lengths,  Figure  19,  indicated  that  for  the  larger 
grid  length,  1AL,  the  solution  oscillates  eastward  of  the  counter 
current  region.  The  smaller  grid  length,  (1/2) al,  had  a  well  defined 
boundary  current  and  the  weaker  counter  current  region.  The  transport 
decreased  rapidly  eastward  of  this  region  with  relatively  little 
oscillation.  Comparing  this  profile  with  the  profile  of  the  modified 
Munk  and  Carrier  solution,  Figure  20,  it  was  observed  that  the  model 
did  in  fact  compare  favorably  with  the  analytic  solution  in  its  pre- 
diction of  the  western  boundary  current  and  the  decrease  in  transport 


39 


CD 

> 

CO 

3: 

>> 

•z. 

_Q 

o 

to 

1 — 1 

to 

h- 

o 

o 

a: 

LU 

a: 

+-> 

»— i 

c 

Q 

O) 

•r— 

X 

C/> 

c 

z. 

fO 

*— 1 

S- 

m 


a. 

p~ 

a 

a> 

»— i 

s- 

a: 

3 

o 

CD 

"T- 
O 

d 


o 


o 


n 


NOIlDNfld  WV3U1S  1VN0ISN3WIQ-N0N 


40 


-I- 
o 

o 


T" 


O 


,_ 

CD 

> 

<X5 

s 

zi 

o 

>> 

1— ( 

XI 

h- 

</> 

O 

to 

LxJ 

o 

or 

a; 

i— i 

Q 

+j 

c 

X 

oj 

•1 — 

zz 

to 

t— 1 

c 

CM 

«T3 

»-    1— 

S- 

zz. 

1— 

t— 1 

o 

Q-  . 

• 

«D 

Q 

i — 

I— i 

C£ 

a> 

C5 

S- 

C7> 


oo 


o  o 

NOIiONflJ  WVBtiJLS  1VN0ISN3WIQ-N0N 


41 


< 

CVJ 

r- 

•o 

, 

c 

fO 

o 

_J 

o 

< 

<o 

S- 

o 

»«- 

c 

o 

•r— 

4-> 

O 

c 

UJ 

3 

s: 

U_ 

t— i 

o*" 

E 
to 

£  -J 

CU 

T   UJ 

%. 

Q 

■*-> 

O 

oo 

s: 

r— 

_i 

A3 

<c 

C 

z: 

o 

o 

•r— 

i— i 

(/» 

oo 

c 

•z. 

cu 

UJ 

E 

s: 

•r— 

t—t 

-a 

a 

1 

O     1 

c 

•O  ZT 

o 

m  o 

z: 

■z. 

cu 

•C 

+-> 

M- 

O 

c 

o 

CO 

•1 — 

i- 

<o 

Q. 

o 

E 

-m 

O 

"p* 

c_> 

cu 

CD 


N0I13NfU  WV3H1S  1VN0ISN3WIQ-N0N 


42 


-a 

c 


o 

4- 

>> 

cr> 

s- 


rC 

o 

•r— 
(/> 

c 

CD 


I 

c 
o 


<u 

4- 
O 

C 

o 
•r- 

S- 

to 
a. 
E 
O 

<_> 


CO 


a; 

s- 

C7> 


A9H3N3  3IJL3NDI  1VNQISN3WIQ-N0N 


43 


a) 

-c 

-t-J 

t/1 

to 

O 

s- 

a 

«a: 

c 

0 

zc 

•r— 

h- 

-!-> 

CD 

a 

2: 

E 

UJ 

3       • 

_l 

U 1 

< 

■z. 

E  co 

1— 1 

<C  >s 

00 

a>  1— 

< 

s- 

CO 

+J  -0 

in  c 

_J 

fO 

<X. 

0) 

^ 

x:  _j 

0 

+j  < 

1—1 

r— 

00 

<+- 

z 

0  s- 

UJ 

0 

2: 

to  M- 

1 — 1 

0) 

0 

r—    s- 

1 

•r-     O) 

2: 

M-    >> 

0 

0    fO 

z 

s 1 

Q_ 

>> 

1—    S- 

fO    ro 

c  -0 

0  c 

•r-     3 

to    0 

C  CQ 

<U 

E    C 

•r-     S- 

-O    0) 

1   +-> 

c  to 

0  <u 

cr> 


s- 

C7> 


NOIlONnj  WV3tflS 


44 


+-> 
co 

CO 

O 

s- 
u 
<c     • 

c 

E    O 

O  t- 
•r-  +-> 
-(->    3 

O  i— 

c  o 
-a 

E  O) 
n3  -r- 
<D  M- 
S-  -r- 

+j  -a 
co  o 


a; 


OJ 


o  s- 
o 
co  ^- 
a> 
r—    s_ 

•i-  a> 
<*-  >> 

O     rt3 
S 1 

i—  i_ 

fO  n3 

C  "O 

O  C 

•i-  3 

CO  O 

C  CO 

CD 
E  c 

•r-  S- 

"O  <D 

I  +-> 

C  CO 

o  a> 


o 
oj 

cu 

s- 

CD 


NOIiONflJ  WV3U1S 


45 


eastward  of  the  counter  current  region.  The  contoured  values  of  the 
stream  function,  Figure  21,  supported  this  observation. 

The  steady  state  solution  generated  by  the  model  could  be  used  to 
acequately  indicate  oceanic  transport.  The  accuracy  of  this  solution 
was  dependent  on  grid  length,  the  more  accurate  solution  being  obtained 
using  the  smaller  grid  length. 
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Figure  21.  Distribution  of  the  Non-dimensional  Stream 
Function  for  1/2aL. 
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V.  CONCLUSIONS 

The  model  was  sufficiently  stable  for  the  lengths  of  integration 
time  used.  The  model  had  to  be  scaled  to  reflect  a  transport  from  zero 
to  100  Sverdrups  with  a  time  step  which  represented  a  one-half  pendulum 
day  to  ensure  stability.  To  reach  steady  state  conditions  it  was  neces- 
sary to  depend  on  this  factor  so  that  the  transients  would  be  damped  out 

The  curl  operation  was  not  accurate  enough  to  produce  the  correct 
results  when  operating  on  the  wind  stress. 

Rossby  waves,  either  topographic  or  inertial,  were  satisfactorily 
described  by  the  model.  The  phase  velocity  was  more  accurately  defined 
with  a  refined  grid.  Although  these  results  were  dependent  on  graphical 
interpretation  this  observed  trend  persisted  through  the  following 
series  of  runs  in  Section  III.  The  errors  were  significant  enough  to 
warrant  future  investigations. 

The  steady  state  solution  generated  by  the  model  compared  favorably 
with  the  modified  Munk  and  Carrier  solution.  To  reach  this  solution  a 
long  integration  time  was  necessary.  The  spin  up  time  was  influenced 
by  mixed  boundary  conditions  which  caused  a  longer  time  to  reach  steady 
state  conditions.  The  transient  solution  the  model  developed  was 
similar  to  that  which  was  reported  by  Gates.  This  dominating  influence 
could  be  more  accurately  described  by  including  some  form  of  bottom 
friction.  These  transients  became  imbedded  in  the  flow  and  were  signi- 
ficant in  establishing  the  counter  current  region  at  steady  state 
conditions.  The  oscillations  observed  in  the  stream  function  profile 
for  larger  grid  length  were  due  to  these  transients  dominating  the 
steady  state  flow. 
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APPENDIX    THE  FINITE  DIFFERENCE  EQUATIONS 

The  specific  details  of  the  derivation  of  the  finite  difference 
equations  that  are  given  below  can  be  found  in  Gait  (1972).  This  deriva- 
tion incompasses  a  combination  of  work  done  by  Winslow  (1967),  Sadourny, 
Arakawa,  and  Mintz  (1967),  and  Williamson  (1969).  By  considering  a 
computational  molecule,  Figure  22,  the  solution  for  each  term  in  equation 
(8)  was  transformed  from  a  surface  integration  to  a  line  integration  by 
using  Stokes  theorem. 


(8) 


(1)  (2)  (3) 

The  value  of  term  (1 )  at  point  L  was  considered  as  an  average  of  the 
potential   vorticity  over  the  molecule  times  the  value  of  the  stream 
function: 
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The  value  of  term  (2)  at  point  L  was  considered  to  be  the  average 
value  over  the  inner  shaded  area  of  Figure  22. 


6v^C  =  6 


L+l  +  ?L-N+L  +  ^L-N  +  ^1  +  ^L-1 


?L+N-1  +  CL+N  _  6CL 


Term  (3)  considered  the  average  values  of  magnitude  and  direction  of  x 
along  each  side  of  the  molecule,  then  stored  as  a  constant  array  since 
it  was  invariant  with  time. 
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The  successive  over-relaxation  procedure  was  used  for  equation   (2) 
where  the  residual   was  computed  by  equation   (9) 

Res     =  vftv  *]  -  %  (9) 

Res  was  determined  by  an  average  value  of  v(  t-  v  ty)     over  the  inner 
hexagon,  Figure  22. 


50 


L+l                    L-N+l                    L-N                   L-l  L+N-l 

Res=4!^-4 ^  +  -r-^ ^  +  -r-4 .  +  ^4 ,-■  +         T 


3 > hL+1  t  hL     hL"N+1  ♦  hL     hL"N  ♦  hL     h1"1  +  hL     hL+fM+  hL 


hL+1+hL/'\hL-N+T+hL/'\hL-N+hL/      Ih1"1   +  hL 


( 


h^1      +hL  \hL  +  N     +     hL 


-  5 


The  amount  of  over-relaxation  was  determined  by  equation  (10). 


*L  =  *L  +  R(ctyL)   R  =  1.48 


(10) 
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